## Figure 3: Effects of Neighborhood Race and Income on Wait Times

## install.packages(c("tidyverse", "fixest"))
## library(tidyverse)
## library(fixest)

## SET WORKING DIRECTORY HERE
## setwd()

## Loading data
## load("dta.RData")

## Figure 3a
## Across-service
race_wait_across <- feols(log_wait_time ~ factor(white_third) |
                            city^open_month^open_year, data = dta, cluster = "geo")
race_wait_across_coefs = as.data.frame(race_wait_across$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Within-service
race_wait_within <- feols(log_wait_time ~ factor(white_third) |
                            city_service^open_month^open_year, data = dta,
                          cluster = "geo")
race_wait_within_coefs = as.data.frame(race_wait_within$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

race_coefs = race_wait_within_coefs %>%
  bind_rows(., race_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(colour = "black"),
        legend.position = "bottom") 

## Figure 3b
## Across-service
inc_wait_across <- feols(log_wait_time ~ factor(inc_third) |
                           city^open_month^open_year, data = dta,
                         cluster = "geo")
inc_wait_across_coefs = as.data.frame(inc_wait_across$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Within-service
inc_wait_within <- feols(log_wait_time ~ factor(inc_third) |
                           city_service^open_month^open_year, data = dta,
                         cluster = "geo")
inc_wait_within_coefs = as.data.frame(inc_wait_within$coeftable) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

inc_coefs = inc_wait_within_coefs %>%
  bind_rows(., inc_wait_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Wait Time \n vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_rect(colour = "black"),
        legend.position = "bottom") 